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ABSTRACT 


Moments of inertia were experimentally determined and longitudinal and 
lateral/directional static and dynamic stability and control derivatives were estimated for a 
fixed wing Unmanned Air Vehicle (UAV). Dynamic responses to various inputs were 
predicted based upon the estimated derivatives. A divergent spiral mode was revealed, 
but no particularly hazardous dynamics were predicted. The aircraft was then 
instrumented with an airspeed indicator, which when combined with the ability to 
determine elevator deflection through trim setting on the flight control transmitter, 
allowed for the determination of the aircraft's neutral point through flight test. The neutral 
point determined experimentally corresponded well to the theoretical neutral point. 
However, further flight testing with improved instrumentation is planned to raise the 
confidence level in the neutral point location. Further flight testing will also include 


dynamic studies in order to refine the estimated stability and control derivatives. 
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I. INTRODUCTION 


A. MISSION NEED 

Unmanned Aerial Vehicles (UAVs) are gaining acceptance as an integral part of the 
operations of today's armed forces. Preceding and during Operation Desert Storm, UAVs 
flew a variety of missions including reconnaissance, targeting for gunfire support, and 
battle damage assessment. Desert Storm provided the first-ever combat test of Unmanned 
Aerial Vehicles by U.S. forces [Ref. 1]. Reviews of lessons learned from Desert Shield 
and Desert Storm reveal the outstanding successes of UAVs. 

There are many examples of effective UAV use during the war. In one instance, a 
commander of a Marine task force was able to monitor UAV imagery of Kuwait as his 
task force approached the city, revealing the exact reaction of the Iraqi forces to Marine 
armor, artillery, and troop movements. The Navy used UAVs to search for mines, spot 
for gunfire support, perform reconnaissance missions for SEAL teams, and search for 
Iraqi Silkworm sites, command and control bunkers, and anti-aircraft artillery sites. The 
Marines quickly reacted to an Iraqi attack into Saudi Arabia observed by a UAV and 
decisively crushed the Iraqi invasion with airborne Cobras and Harriers. The Army 
provided their Apache pilots with route reconnaissance acquired from UAVs shortly 


before the Apache missions. Spotting for air strikes and naval gunfire support 


became so successful that Iraqi soldiers were seen attempting to surrender to UAVs as 
they flew overhead [Ref. 2]. 

UAVs have many advantages over manned aircraft which help account for their 
effectiveness. First, the cost of a UAV is a very small fraction of the cost of a manned 
aircraft. The Pioneer, for example, costs approximately $500,000 for the aircraft and 
$500,000 for the onboard camera, bringing the total cost of the package to a mere one to 
two percent of the cost of most manned tactical aircraft. UAVs are also extremely flexible 
with respect to launch platform. They can be launched from small fields, truck beds, and 
practically any ship in the Navy inventory. UAVs are frequently very hard to detect with 
radar or infrared systems due to their small size, composite construction, small engines 
(sometimes electric motors), and their slow speed. UAVs also have an advantage from 
being unmanned. The aircraft is not limited by the "g" tolerance of a pilot or by pilot 
fatigue. Finally, the best advantage of all is that when a UAV crashes or 1s lost to enemy 
fire, there is no search and rescue mission required, no prisoners of war taken, and no loss 
of life. 

UAVs are not without their problems, however. Acquisition and support programs 
are relatively new and underdeveloped. This results ina UAV force that is relatively small 
in number of aircraft, and small and inexperienced in terms of personnel. The Pioneer 
showed signs of these underlying problems during Desert Storm. Six Pioneers were 
damaged badly enough to require return to the factory for repair due to no intermediate 


level maintenance facilities being available. Five Pioneers were lost due to mechanical 


malfunction or operator error (Ref. 2]. Due to the small number of available aircraft, 
losses are very costly and operator errors need to be avoided if possible. Thorough and 
frequent training through simulation can keep operators performing at their peak. 
B. STATEMENT OF OBJECTIVE 

To provide realistic simulation for training, accurate aerodynamic data for the 
aircraft must be available. Aerodynamic parameters of the aircraft can be determined from 
its physical characteristics, wind tunnel tests, and flight tests of actual or scale models. It 
Is the objective of this work to provide the ground work for determining whether flight 
test can effectively be used to accurately estimate the static and dynamic stability and 
control characteristics of a fixed wing UAV. Flying-qualities parameters were estimated 
using analytic techniques, and the resultant flight dynamics were simulated to provide 


expected behavior for future test flights. 


Il. BACKGROUND 


A. EVOLUTION OF AN AIRCRAFT 

For any aircraft to get from an idea to an actual flying vehicle, properties such as 
stability and control derivatives, which determine the flying characteristics of the aircraft, 
must be determined. These derivatives are used to size control surfaces, design flight 
control systems, and program training devices such as simulators. There are typically 
three ways to determine or estimate derivatives. 

The first method of determining an aircraft's derivatives begins early and continues 
throughout the design process. This step involves determining derivatives mathematically 
from physical characteristics of the aircraft. Requiring little more than a few 
well-educated engineers and some calculators or desktop computers, this method is 
relatively simple. Derivatives can be reasonably approximated, but must be refined 
through other methods. 

The second method of determining derivatives is through aerodynamic 
measurements of wind tunnel testing. This method is more complicated than simple pencil 
and paper calculations due to the necessity for model making and wind tunnel operation, 
but much better results can be achieved. Derivatives must still be refined, however, due to 


factors such as scale effects and interference from wind tunnel walls and 


supporting hardware. Also, dynamic effects are often difficult to properly account for in 
most wind tunnels. 

The final step approach to determining the derivatives of an aircraft is through flight 
test. This method is by far the most expensive and complicated way to determine the 
aircraft's derivatives, but it is also the most precise and complete. Scaled flight test 
provides a viable option for this third method. 

B. PIONEER UAV 

In June 1982, Israeli forces very successfully used UAVs as a key element in their 
attack on Syria. Scout and Mastiff UAVs were used to locate and classify SAM and AAA 
weaponry and to act as decoys for other aircraft. This action resulted in heavy Syrian 
losses and minimal Israeli losses. A year and a half later, the U.S. Navy launched strikes 
against Syrian forces in the same area with losses much higher for the Navy than those of 
the Israelis [Ref. 3]. 

The Commandant of the Marine Corps, General P. X. Kelly, recognized the 
effectiveness of the Israeli UAVs. Secretary of the Navy John Lehman then initiated 
development of a UAV program for the U.S. Anxious to get UAVs to the fleet, Secretary 
Lehman stipulated that UAV technology would be off-the-shelf [Ref. 4]. After the 
contract award to AAI Corporation of Baltimore, Maryland for the Pioneer UAV, 
developmental and operational testing took place concurrently. This approach resulted in 


quick integration of the Pioneer into the fleet. Unfortunately, such quick integration into 


the fleet can result in problems identified during operational use which had not been fully 
explored in the test and evaluation process. 

The UAV Office at the Pacific Missile Test Center (PMTC, now the Naval Air 
Warfare Center, NAWC, Weapons Division, Pt. Mugu) was tasked with Developmental 
Test and Evaluation of the Pioneer. Testing revealed the following concerns which 


warranted further investigation [Refs. 5,6]: 


l. discrepancies in predicted with flight-tested rate of climb, time to climb, 
and fuel flow at altitude; 


2. | apparent autopilot-related pitch instability; 
3. tail boom structural failure; 
4. | severely limited lateral control; 


5. slow pitch response causing degraded maneuverability at high gross 
weights; 


6. insufficient testing to determine the effects of the new wing on flight 
endurance. 


The Target Simulation Laboratory at Pt. Mugu was tasked to develop a computer 
simulation of the Pioneer in order to provide cost-effective training for pilots. 
Aerodynamic data were needed to provide the stability and control derivatives necessary 
for the simulation as well as to answer questions concerning basic flying qualities of the 
Pioneer. 

In order to provide support to the research being done at Pt. Mugu and to provide 
for future UAV project support, a research program was begun at the Naval Postgraduate 


School (NPS). An instrumented half-scale radio-controlled model of the Pioneer was used 


for the research at NPS. Research performed included wind tunnel tests, flight tests, and 
numerical modeling. 

Initial NPS research on the Pioneer, performed by Capt. Daniel Lyons, involved a 
computer analysis of the Pioneer in its original configuration and with a proposed larger 
tail. A low order panel method (PMARC) was used for the aerodynamic analysis. Static 
longitudinal and directional stability derivatives, the neutral point, and crosswind 
limitations were calculated. Drag polars were constructed using the component buildup 
method for profile drag, and drag reduction measures were considered [Ref. 7]. 

In conjunction with Capt. Lyons work, Lt. James Tanner conducted wind tunnel 
tests to determine propeller efficiencies and thrust coefficients for drag studies [Ref. 8]. 
Lt. Tanner also conducted flight tests to determine power required curves and drag polars 
(Ref. 8]. Capt. Robert Bray later conducted wind tunnel tests of a 0.4-scale model at 
Wichita State University to determine static stability and control derivatives [Ref. 9]. 
Aerodynamic data obtained by Capt. Lyons and Bray have been supplied to PMTC to be 
used for simulation. 

Lt. Jim Salmons performed initial flying qualities flight testing rens an onboard data 
recording system in order to determine static stability parameters. Unfortunately, 
vibration problems with the onboard recorder rendered much of the data unusable (Ref. 
10]. 

Following up on Lt. Salmons' work, Lt. Kent Aitcheson installed the CHOW-1G 


telemetry system, designed by Lt. Kevin Wilhelm, in an attempt to alleviate the vibration 


problem experienced by Lt. Salmons. The new flight test configuration was used to test 
static longitudinal and lateral-directional stability characteristics of the Pioneer. The 
vibration problem experienced by Lt. Salmons was overcome, though not enough data 
were acquired for a complete and thorough analysis of the Pioneer's characteristics. Much 
insight was gained, however, concerning instrumentation. Resolution needed to be 
improved for flight control position indication [Refs. 11,12]. 

Lt. Paul Koch conducted further flight tests of the Pioneer with the CHOW-1G 
telemetry system. Static longitudinal stability results from the flight tests correlated well 
with theoretical predictions and with simulations of a full-scale Pioneer. Electromagnetic 
interference with the flight control system at the test site resulted in loss of the half-scale 
model Pioneer before further data collection and analysis could be performed [Ref. 13]. 

C. PARAMETER ESTIMATION 

Parameter estimation is used to derive stability and control derivatives from dynamic 
flight test data. Lcdr. Robert Graham successfully used the Modified Maximum 
Likelihood Estimation (MMLE) technique with MATLAB software to analyze simulated 
and actual flight test data of several aircraft. Analysis of simulated UAV data revealed the 
effect of signal-to-noise ratio on the estimator, and the need for proper control-surface 
excitation for a particular response [Ref. 14]. Cdr. Patrick Quinn analyzed flight test data 
from the Marine Corps BQM-147 UAV using both MMLE and a more robust non-linear 
model, pEst, and compared the results of the two approaches. Noise was found to be a 


problem when the system response was in the same general frequency as the noise. 


Limited available data and noise at the frequency expected for the system's response 
prevented a successful resolution of all stability and control derivatives of interest. The 
pEst model was found to be the estimator of choice when aircraft maneuvers exceed what 
is generally considered reasonable for linear approximation of flight dynamics [Ref. 15]. 
While previous work with UAVs at NPS has at times been frustrating, much has 
been learned, especially concerning the most challenging method of aircraft analysis, flight 
test. This work initiates the implementation of the lessons learned from previous work 
into dynamic simulation and flight testing of a generic UAV in order to properly prepare 
the UAV lab at NPS for further support of future projects and to demonstrate the value of 


scaled UAV flight test to the fleet. 


Il. THE AIRCRAFT 


A. GENERAL DESCRIPTION 

The "Bluebird", shown in Figures 1-3, is a high-wing tricycle-gear radio-controlled 
airplane. It is constructed of wood, foam, composites, and metal. It is powered by a 
Sachs-Dolmar 3.7-cubic-inch two-stroke gasoline engine which drives a 24 inch 
two-bladed wood propeller. It is controlled by a nine-channel pulse-code-modulated 
Futaba radio operating at 72.710 MHz. To enhance reliability, the Bluebird has two 
receivers which share control of the aircraft. The left receiver controls the left aileron, 
elevator, and flap, and engine ignition and onboard electronics package cut-off, while the 
right receiver controls the right aileron, elevator, and flap, and rudder, nose-wheel 
steering, and throttle. The Bluebird can fly within visual range for approximately 1.5 
hours. Table 1 describes physical specifications of the Bluebird. 
B. STABILITY DERIVATIVES 

The initial estimates of the stability derivatives of the Bluebird were made using the 
physical characteristics of the aircraft such as airfoil data , geometric measurements, 
relative positions of aircraft components, mass, and weight (Refs. 16-19]. The assumed 
flight condition with the associated aircraft configuration is described in Table 2. 


Nondimensional stability and control derivatives estimated are shown in Table 3, and 
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dimensional stability and control derivatives estimated are shown in Table 4. MATLAB 


programs written for the physical and derivative calculations appear in Appendix A. 


TABLE 1 
SPECIFICATIONS 

E ee e n n ns 9.84 ft 
ANO —-— s E 3.04 ft. 
Wing Airfoil (est.) .......................... GO 769 
Horizontal Stab. Airfoil (est.) ........... NACA 4412 
SMS) A e erede 22.38 ft.? 

. EE. H. H 7 wa e's 4.701 ft.* 
SE ES oL Lee rem eR du» T fs 
CEN a ooo nig eds 1.802 ft 
Oc oux m. ND sco As ss te tl 1.281 ft 
HESS 2 OEC EN act a 12.42 ft 
DM AE LL THT E 3.67 ft 
DENS A Ri Pe cco AM esee 1.21 ft 
E 6.89 
PE H a NEL is sse 2.86 
SENE a bs 1.14 

A ca 0.6274 
A cesis eee 0.0247 
E a O es 0.0023 
HE a 5.381 ft. 
E ti 5.381 ft. 

TABLE 2 

FLIGHT CONDITION/AIRCRAFT CONFIGURATION 
ANE ITem a a a D 57.79 Ibs. 
O o o Es 12.58 slug'ft' 
IEEE OE cc gnal 13.21 slug*ft? 
Eas oes v 19.99 slug*ft' 
A a e 60 mph/88 ft/sec 
AM o nel E eye. sss.» aga - 800 ft MSL 
BERS a E 0.002327 slugs/ft? 
Center of Gravity ...................... 27% of m.a.c. 
O E a E E 0.2866 
A a ee ee 3.8 degrees 
SOYER ee a. 1.6 degrees 





Figure 1 
Side View 


Figure 2 
Top View 





Figure 3 
Front View 


TABLE 3 
NONDIMENSIONAL DERIVATIVES 


BEES V 0.2866 OT C 0.0358 
E TT 4.1417 Um S 0.0311 
E ocu 1.5787 A ies 0.1370 
Zo nu -1.0636 S EEE -4.6790 
— 3.9173 Ce -11.6918 
AAA suus 0.4130 —— 0.0650 
B -1.2242 A -0.3100 
CA 0.0484 Ce A -0.0330 
Em, om 0.0000 Gm a -0.0358 
CI ERN -0.3579 COMME ess 0.0967 
M m m oA -0.0526 C; mmm. 0.0755 
EHE 0.0000 CR i -0.0258 
A 0.2652 A mus 0.0697 
"BA -0.0326 E 0.0028 
Soa 0.0000 (Cm m... 0.0000 
TABLE 4 


DIMENSIONAL DERIVATIVES 


EW W......... -0.0914 /sec P cm 16.7894 fUsec 
Pie L, -7.2961 ft/sec? V RE -0.7312 /sec 
EM ..... -468.9852 ft/sec? LA O 1.8146 ft/sec 
ZA 4.5027 ft/sec WE uuu. -46.368 ft/sec? 
EN uuu 0.0000 /ft*sec MO rok. -29.2559 /sec* 
PP -1.3178 /sec ee -3.2928 /sec 
ES ees -33.6730 /sec* YOUR E -34.8021 ft/sec? 
ae 0.0000 ft/sec MORE s. 0.7663 ft/sec 
IE... 0.0000 ft/sec* — -7.8282 ft/sec? 
_ RE -6.5787 Isec' Deren -5.0281 /sec 
pnm... 1.0613 /sec Dome 52.7966 /sec? 
ees 0.5589 /sec? im es 6.0593 /sec? 
BENE O -0.3167 /sec D eto a -0.4647 /sec 
o RE -3.2375 /sec? ce -4.0900 /sec? 
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C. MOMENTS OF INERTIA 

Having accurate moments of inertia is critical to ensuring the accurate prediction of 
aircraft dynamics. Direct calculation of a model's moments of inertia by consideration of 
the contributions made by individual parts is impractical and inaccurate. Determination of 
the moments of inertia by test is much more practical and precise. The changes in a 
model's moments of inertia due to addition or subtraction of equipment or structure can be 
calculated directly, thereafter. 

In the determination of moments of inertia by test, the aircraft is hung from the 
ceiling and swung. Using the period exhibited and the principles of compound pendulums, 
the moments of inertia of the model can be extracted [Ref. 20-21]. In order to calculate 
the moment of inertia about all axes, the model must be hung from the ceiling and swung 
three different ways, each such that as the aircraft swings, it is rotating about the axis of 
interest. The Bluebird was hung by chain and swung as pictured in Figures 4-6. 
Specifications for the geometry of each test can be found in Appendix B. 

Reference 21 provides equation (1) for calculating the moment of inertia of a 


swinging model. 





- 5 32 
_ Wms+sZmMsS 2 WsZs 2 WuZm 
[= 4n2 Pm+s 7 au Ps — Tg (1) 


W is the weight, Z is the distance from pivot to center of gravity, p is the period, and g is 


the gravitational constant. M and $ are subscripts designating either the model or the 


support. 
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Figure 4 
I. Test (not to scale) 





Figure 5 
I, Test (not to scale) 





Figure 6 
[,, Test (not to scale) 


It was determined that swinging the support (the chains) in the configuration it 
would be in when supporting the model would not be possible, since the chains would not 
maintain their positions without the model in place. Equation (1) was therefore 
manipulated in order to treat the chains as long slender rods and to calculate their 
moments of inertia as such [Ref. 22]. 

The new form of equation (1) is 

LE E PEE Q) 
where Ls is the length of a chain and the summation is taken over all chains (four in this 


case). In particular, there were two "iong" chains and two "short" chains, all having a 


weight per unit length o. Appropriate substitutions were made to yield 


= =? 
_ [Wu +20(Lshorr + Lion) mes Hu ZM 20 ([ y 2 2 3 
is TY E AA MST g 237 L snort + Lion (3) 


Having the equation in this form fixes the values for all variables except 
Žas, Zu, and Puss. These three variables were measured for each of the three 
configurations and calculations were made. Four periods were timed during the swing 
tests, and the tests were repeated ten times. The moments of inertia calculated are shown 


in Table 2. 
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IV. PREDICTED DYNAMICS 


A. GENERAL 

When flight testing an aircraft, it is important to attempt to predict the results of the 
testing prior to the actual flights. The information provided by the predictions can be used 
in briefing the pilot as to what to expect from the aircraft and also to avoid any potentially 
dangerous flight regimes. Longitudinal and lateral-directional dynamics of the Bluebird 
have been predicted using the stability and control derivatives estimated in chapter three 
along with computational methods based upon the six equations of motion as described in 
references 18 and 23. Computational programs were written in MATLAB and appear in 
Appendix C. 
B. LONGITUDINAL DYNAMICS 

The MATLAB program named "longnat.m" uses the full (4x4) longitudinal plant 
and the MATLAB "impulse" function to determine the short and long period natural 
responses. A diary of the values computed and output by the program provides 
eigenvalues, damping ratios, and damped and undamped natural frequencies for both 


modes as shown in Table 5. 
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TABLE 5 


LONGITUDINAL DATA 
Short Period Long Period 
Eigenvalues -5.083 +/- 4.8611 -0.037 +/- 0.400i 
Damping ratio 0.723 0.093 
Undamped Natural Frequency 7.03 rad/sec 0.401 rad/sec 
Damped Natural Frequency 4.86 rad/sec 0.399 rad/sec 


Figures 7 and 8 show the combined short and long period natural response. It can be seen 
that the short period response is almost completely damped out after only two seconds. 
Response beyond two seconds is primarily long period. The phugoid mode is seen to be 
lightly damped. 

The short and long period response to a unit step elevator input is determined by 
using the full (4x4) longitudinal plant and the MATLAB "step" function in the MATLAB 
program named "stepper.m". Figures 9 and 10 show the short and long period response 
to a step input. In the first plot, the long period can be seen to be much more heavily 
excited than the short period mode. Again, it can be seen that the short period response is 
almost completely damped out after about two seconds. Due to held elevator input, the 
long-term response is to trim to a new angle of attack, while pitch rate dies out. 

The MATLAB program named "n step.m" uses the full (4x4) longitudinal plant and 
the MATLAB "step" function to determine the normal acceleration or load factor 
response to a unit step elevator input. Figure 11 shows the short period normal 
acceleration response to a unit step input. It should be noted that normal acceleration is 


defined opposite in sign to load factor. The long-period response is seen to be much 
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Figure 7 
Natural Response (short time) 
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Figure 8 
Natural Response (long time) 
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Figure 9 
Response to Step Input (short time) 
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Figure 10 
Response to Step Input (long time) 
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more critical than the short-period response in generating large load factors for this 


aircraft. 
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Figure 11 
Normal Acceleration Response to Step Input 


The longitudinal response to given initial conditions is determined using the full 
(4x4) longitudinal plant and the MATLAB "initial" function in the MATLAB program 
named "homogen.m". The initial conditions ( u/U = .34, alpha = 5 deg, q = 8.8 deg/sec, 
theta = -0.8 deg) are provided for all states from the step input results after initial 
dynamics have died out (approximately 15 seconds). Figures 12 and 13 show the 
longitudinal response to the initial conditions. The response in angle of attack is small 
while a significant pitch rate is developed in both the short-period and long-period 
responses. The long-period is again seen to be lightly damped. 

The MATLAB program named "doublet.m" uses the full (4x4) longitudinal plant 


and the transition matrix (using the MATLAB "exponential" function) to find the 
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Figure 12 
Homogeneous Response (short time) 





Figure 13 
Homogeneous Response (long time) 
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response to a unit pitch doublet input at approximately the short period damped natural 
frequency. Figure 14 shows the input pitch doublet, and Figure 15 shows the response. 
Though the doublet is commonly used to observe only the short-period response, the 
angle of attack response indicates the long-period response is also excited. 

The response (transfer function gain and phase) to a harmonic elevator input is 
found using the full (4x4) longitudinal plant and the MATLAB "bode" function in the 
MATLAB program named "sp bode.m". Figures 16 and 17 show the gains and phase 
shifts from an elevator frequency sweep. The angle of attack gain can be seen decreasing 
from a peak at a very low excitation frequency; the short-period response is masked by the 
long-period gain. In actuality, the low-frequency response is of little interest. Pitch rate 
shows a maximum gain at approximately 2.4 radians/second corresponding to an elevator 
cycle period of approximately 2.6 seconds. 

The MATLAB program named "n bode.m" uses the full (4x4) longitudinal plant 
and the MATLAB "bode" function to find the normal acceleration or load factor response 
(transfer function gain and phase) to a harmonic input. Figures 18 and 19 show the 
normal acceleration gains and phase shifts from an elevator frequency sweep. The short 
and long-period damped natural frequencies can be observed at 4.86 radians/second and 
0.399 radians/second respectively where their respective gains peak. While the gain 
demonstrated at the long-period damped natural frequency is quite significant, the aircraft 
would not be expected to be excited at that frequency. An excitation near the 


short-period damped natural frequency is much more likely. Flight test pilots and 
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Figure 14 
Unit Elevator Doublet Input 
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Figure 15 
Doublet Response 
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Figure 16 
Frequency Response (Gain) 
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Figure 17 
Frequency Response (Phase) 
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Figure 18 
Normal Acceleration Response to a Harmonic Input (Gain) 
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Figure 19 
Normal Acceleration Response to a Harmonic Input (Phase) 
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engineers should be aware of the fact that excitation at that frequency will produce 
approximately 0.15 "g's" per degree of elevator deflection or that it takes a harmonic 
amplitude of about 6.7 degrees of elevator to produce one "g" of acceleration. The 
Bluebird has a maximum of 15 degrees of up elevator and 12 degrees of down elevator 
available. This would equate to approximately -0.8 to *3.2 "g's", which is easily tolerated 
structurally. 
C. LATERAL-DIRECTIONAL DYNAMICS 

Undamped natural frequency, damped natural frequency, damped natural period, 
and damping ratio for the Dutch roll mode; time constant for the roll mode; and time 
constant and time to double or half for the spiral mode are calculated by the MATLAB 


program named "lat dir.m" using the full (4x4) lateral-directional plant. Values calculated 


are described as follows. 


Dutch roll mode: 


Undamped natural frequency 2.65 rad/sec 

Damped natural frequency 2.62 rad/sec 

Damped natural period 2.40 sec 

Damping ratio 0.148 
Roll mode: 

Time constant 0.195 sec 


Spiral mode: 
Time constant -29.28 sec 
Time to double amplitude 20.29 sec 


The Dutch roll mode is observed to be lightly damped. Also, the spiral mode is divergent 


with a fairly short time to double bank angle. 
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The MATLAB program named "rudkick.m" uses the full (4x4) lateral-directional 
plant and the MATLAB "impulse" function to find the response to a unit rudder impulse. 
The program also finds the bank angle at the end of 100 seconds which is approximately 
four degrees in the direction of the rudder kick. Figure 20 shows the unit rudder kick 
response in bank and sideslip angles. Sideslip lags bank angle by approximately 0.6 
seconds or 80 degrees of phase and is approximately 60 percent larger in magnitude. 

The response to an initial sideslip is determined by the MATLAB program named 
"sideslip.m" using the full (4x4) lateral-directional plant and the MATLAB "initial" 
function. Initial conditions for sideslip and bank angle are provided by a steady-sideslip 
condition where bank angle is ten degrees, and sideslip angle is 13.7 degrees. Figure 21 
shows the homogeneous response to the initial sideslip. The lightly damped Dutch roll 
mode can be seen by observing the oscillations of both bank angle and sideslip angle, while 
the divergent spiral mode can be observed by the increasing bank angle. 

The MATLAB program named "roll.m" uses the full (4x4) lateral-directional plant 
and the MATLAB "bode" function to find the roll angle response (transfer function gain 
and phase) due to a harmonic aileron input. Figures 22 and 23 show the roll angle gains 
and phase shifts from the harmonic aileron input. The maximum high-frequency roll angle 
gain is seen to occur at approximately 2.8 radians/second which corresponds to an aileron 


input period of approximately 2.2 seconds. 
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Figure 20 
Unit Rudder Kick Response 
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Figure 21 
Homogeneous Response to an Initial Sideslip 
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Figure 22 
Roll Angle Gain due to a Harmonic Aileron Input 
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Figure 23 
Roll Angle Phase Shift due to a Harmonic Aileron Input 
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D. REMARKS 

Analysis has been conducted for the most common modes. The Bluebird 
demonstrates no particularly dangerous flying qualities. Particular characteristics of the 
Bluebird's behavior worth noting include a very heavily damped short period longitudinal 
mode, a lightly damped phugoid mode, a lightly damped Dutch roll mode, and a divergent 
spiral mode. Additional analysis could be easily done by analogy to those modes all ready 
analyzed. Further analysis might include, for example, the response to a step aileron 


input, aileron doublet, or elevator impulse ("stick rap"). 
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V. FLIGHT TEST 


A. TEST DESCRIPTION 

The relative positions of an aircraft's center of gravity and neutral point are critical 
to the aircraft's longitudinal stability and handling qualities. In order for a conventional 
(aft tail) aircraft to be statically stable, its center of gravity must be forward of its neutral 
point. The farther the center of gravity is in front of the neutral point (relative to the 
chord length), the more statically stable the aircraft is. As the center of gravity is moved 
aft beyond the neutral point, the aircraft becomes statically unstable and more 
maneuverable. 

As an aircraft's center of gravity is moved aft toward the neutral point, less and less 
change in elevator trim is required to achieve steady, level flight for a given change in 
airspeed. This fact can be used to experimentally determine the aircraft's neutral point. 
Reference 23 provides equation (4) which describes the relationship between change in 
pitching moment coefficient with change in coefficient of lift and the required change in 
elevator deflection with change in coefficient of lift. 

== VAS (4) 


As the center of gravity moves aft and approaches the neutral point, dC, /dC, approaches 


zero. Since V, and C, are constants, dóe/dC, must also approach zero. 
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In order to determine the neutral point by flight test, the aircraft must be flown at 
various centers of gravity. At each particular center of gravity, the aircraft is trimmed at 
several airspeeds (coefficients of lift). Each airspeed and the corresponding elevator 
deflection are recorded. Plots are then generated showing elevator deflection versus 
coefficient of lift at each center of gravity tested. A line is then best fit through the data 
points for each center of gravity. The slope of this line represents dóe/dC, for that 
particular center of gravity. Finally, dóe/dC, versus center of gravity is plotted. A best fit 
line is then drawn through these data points. Using the best fit line just drawn, the center 
of gravity where dóe/dC, equals zero is the aircraft's neutral point. 

B. INSTRUMENTATION 

In order to acquire the airspeed of the Bluebird, a simple, commercially-available 
airspeed indicator was installed. The Digicon TT-01 Tele Tachometer/ASI senses the 
spinning of a small wind-driven propeller blade using a cadmium disulphide optical sensor. 
The frequency of the changes in light intensity sensed is transmitted to a hand-held 
receiver which converts the frequency to airspeed which can be read directly in real-time. 
The manufacturer's claimed accuracy is +/- 0.5 feet per second. The airspeed indicator 
was installed on the left wing of the Bluebird by mounting it on the end of a boom which 
extended approximately 18 inches (approximately 80 percent of c-bar) in front of the 
leading edge of the wing. 

Measuring elevator trim was done in an indirect manner. The flight control 


transmitter has a digital display which can show elevator trim on a generic scale of -100 to 
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100. Prior to flight test, elevator deflection in degrees was measured for various trim 
settings, and the relationship between transmitter displayed trim number and actual 
degrees of elevator deflection was determined. This calibration allowed for the 
determination of elevator deflection in degrees during post flight analysis. 

C. TEST PROCEDURES 

In order to fly the Bluebird at various centers of gravity, an access panel in the aft 
fuselage of the aircraft was modified to accept added weights. The aircraft was then 
configured with various amounts of weight and the location of the center of gravity 
determined for each configuration. Centers of gravity were determined by weighing each 
wheel with a strain-gage balance. 

A total of seven flights were flown at various centers of gravity. Flights were kept 
short in order to minimize the shift in the center of gravity due to fuel burn. The location 
of the center of gravity was determined both before and after the flight, and the average 
center of gravity was used for calculations. Shifts in center of gravity during testing were 
kept to one percent or less of the wing chord. 

Each flight consisted of 12 passes at various airspeeds. Airspeed and trim setting 
for each pass were recorded for post flight analysis. Additionally, air pressure, air 
temperature, and aircraft weight were noted in order to calculate coefficient of lift. 

Post flight analysis of the data was conducted in accordance with the procedures 
described in Section A, Test Description. The resulting graphs are shown in Figures 


24-31. 
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Figure 24 
Flight 1 of 7 
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Figure 25 
Flight 2 of 7 


37 


C.G. = .2868 
04 MAR 94, Flight 1 


rees of Elevator 





Slope = -7.461591 





Figure 26 
Flight 3 of 7 
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Figure 27 
Flight 4 of 7 
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Figure 28 
Flight 5 of 7 
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Figure 29 
Flight 6 of 7 
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Figure 30 
Flight 7 of 7 
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Figure 31 
All Flights 
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D. RESULTS 

The neutral point estimated by flight test was found to be about 54 percent of the 
mean aerodynamic chord. The neutral point estimated through conventional calculations 
was approximately 52.8 percent of the mean aerodynamic chord. 

While such a close match between the neutral point locations found by each method 
is highly encouraging, one must not overlook the scatter in the data presented in Figure 
31. Assuming a normal distribution of the data, there is a 68 percent confidence that the 
experimentally determined neutral point is within five percent of the wing mean 
aerodynamic chord of the location determined experimentally. Such data scatter is 
unacceptable for accuracies required, and the tests will be repeated when the improved 
instrumentation under development comes on line. 

For the flight testing done, two sources of potential error are of particular interest. 
First, it was very difficult to ensure a perfectly level pass of the aircraft at each particular 
airspeed. The pilot was positioned on the ground and the plane was flying approximately 
one hundred feet above him. This is not an optimum vantage point for the pilot. Ideally, 
the pilot would like to be elevated (such as in a tower) such that the aircraft is at eye level 
for each pass. Future plans include using an altitude-hold autopilot to maintain the desired 
trim conditions. Also, due to wind gusts and possibly small unintentional changes in 
aircraft attitude, the airspeed was not always steady, and therefore somewhat difficult to 
read consistently. This second problem will be overcome by the improved data acquisition 


system. 
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VI SUMMARY 


A determination of the moments of inertia by the component contribution method 
was decided to be too cumbersome for analysis of the complete aircraft. Moments of 
inertia were found through a compound pendulum analysis, and appear reasonable. Future 
changes to the configuration of the aircraft can be accounted for by considering the 
contribution of the component added, removed, or relocated. 

Initial estimates of stability and control derivatives were made by conventional 
aircraft-design-type methods. Such methods involve considering characteristics of the 
aircraft such as lift-curve slopes of airfoils and the locations of particular aircraft 
components relative to each other. Programs were written in MATLAB in such a manner 
that future changes to the aircraft or different flight conditions can be accounted for 
quickly and accurately. The initial estimates of the derivatives can be used in the 
preliminary design of a flight controller. 

The initial estimates of the stability and control derivatives were used to predict the 
dynamic response of the aircraft to several different excitations. MATLAB programs 
were written to predict the dynamics, and any future modifications to the aircraft or flight 
conditions can be accounted for easily. Several characteristics of the aircraft's responses 
are worth noting. First, the short period longitudinal response was found to be heavily 


damped. Little or no evidence of the short period response can be observed beyond 
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approximately two seconds. Second, the long period longitudinal response was found to 
be lightly damped. The long period response can still be observed relatively easily after 
100 seconds. Similarly, the Dutch roll mode is also fairly lightly damped and can be 
observed easily for ten to 15 seconds. Finally, the spiral mode was found to be divergent 
with a time to double bank angle of approximately 20 seconds. While this mode is not 
particularly dangerous for a radio controlled aircraft which is flown visually (open loop), 
the pilot should be aware of this tendency of the aircraft in order to avoid trim problems. 
It is also recommended that flight controller designs incorporate bank angle feedback for 
wing leveling. 

The neutral point of the aircraft was found to be 53 to 54 percent of the wing mean 
aerodynamic chord by two different methods. The data collected through flight test were 
somewhat scattered; it is estimated that the neutral point location is known within +/- five 
percent of the wing mean aerodynamic chord with a confidence of about 68 percent. 
Further flight testing with improved instrumentation is recommended to raise the 
confidence of the neutral point estimation. It is recommended that flight tests with 
improved instrumentation continue to verify or update all predicted stability and control 


derivatives. 
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APPENDIX A: MATLAB PROGRAMS USED TO ESTIMATE 
STABILITY AND CONTROL DERIVATIVES 


1. blbrdfcl.m 


% Eric J. Watkiss 

% AA0810 Thesis 

% File for Bluebird data which change with flight condition 
% blbrdfcl.m 

% Last Update: 02 MAR 94 


g = 32.174; % Accelleration due to gravity 
Wlmg - 24.02 % Weight on left main in lbs 
Wrmg = 23.91 % Weight on night main in lbs 
Wng = 11.07 % Weight on nose gear in lbs 
Umph = 60 % Flight speed in miles per hour 
Ufps = 88.0 % Flight speed in feet per second 
rho = .002327 % Air density in slugs/(cubic ft) 
Ixx = 12.58 % Moment of inertia about x-axis 
Iyy = 13.21 % Moment of inertia about y-axis 
Izz = 19.99 % Moment of inertia about z-axis 
Ixz = 0 9o Assumed!!! 1H TH TH HE HE IT 
LD - 8; % Lift to drag ratio 

thetanaut = 0; % Initial pitch angle 
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2. bluebird.m 


% Eric J. Watkiss 
95 AA0810 Thesis 


95 File for Bluebird data which are fixed 


9/5 bluebird.m 


% Last Update: 02 MAR 94 


ac = .479; 

al = 3; 

% 

alphaOl = -6.5*p1/180; 
% 


% 

b = 12.42; 

bt = 3.67; 

bv = 1.208; 
cbar = 1.802; 
% 


CLalphaafw = 5.443; 
% 

% 

CLalphaaft = 5.587; 
% 

% 

CLalphaafv = 2* pi; 
% 

% 

CMac = -.06; 

% 

ct = 1.281; 

c4tail = 7.878; 

% 

% 

c4wing = 2.497; 

% 

da0dde = 625; 

% 

% 
daOddr 
% 

% 
deda = .4; 


675; 


90 Aileron chord in ft. 

% distance from centerline to 
inner edge of aileron in ft. 

% a.o.a. for zero lift (radians)ao = 6; 
distance from centerline to 
outer edge of aileron in ft. 

% Span of wing in ft 

% Span of horizontal tail in ft. 

% Height of vertical tail in ft. 

% Mean aerodynamic chord (m.a.c.) 
in feet 

% Lift curve slope of wing 
airfoil (GO 769) in per 
radian 

% Lift curve slope of horizontal 
tail airfoil (NACA 4412) in per 
radian 

% Lift curve slope of vertical 
tail airfoil (flat plate) in per 
radian 

% Coefficient of moment about 
aero. ctr. (GO 769) 

% m.a.c. of horizontal tail in ft. 
% Location of quarter chord of 
horizontal tail in feet from 

firewall 

% Location of quarter chord of 
wing in feet from firewall 

% Section flap effectiveness 
for 33% flap (elevator) 
Abbott and Doenhoff p. 190 

% Section flap effectiveness 
for 38% flap (rudder) 
Abbott and Doenhoff p. 190 

% Downwash angle derivative 
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% 

Df = 1; 

e0 = Q; 

ee =.8; 

g = 32.174; 
hac = .245; 
% 


it = 4.83 *p1/180; 
lewing = 2.047; 
% 

letail = 7.557: 
% 

% 

mg = 37.595/12; 
% 


ng = .75/12; 
% 

S = 22.380; 
Sr = .547; 
St = 4.701; 
Sv = 1.277; 
Wf = .67; 
ybar = b/4; 
zv =; 

% 
Zwt=.S; 
% 


estimated from Perkins/Hage 

% Depth of fuse. in ft. 

% Assumed epsilon naught 

% Assumed span efficiency factor 

% gravitational constant 

% Location in percent chord of 
aero. ctr. (NACA 4412) 

% Incidence angle of hor. tail 

% Location of leading edge of wing 
in feet from firewall 

% Location of leading edge of 
horizontal tail in feet from 
firewall 

% Location of main gear in ft 
from firewall 

% Location of nose gear in ft 
from firewall 

% Reference (wing) area in sq. ft. 

% Rudder area in sq. ft. 

% Horizontal tail area in sq. ft. 

% Vertical tail area in sq. ft. 

% Width of fuse. in ft. 

% Spanwise location of m.a.c. 

% Vert. tail height to m.a.c. 
(estimated) 

% Verticle height of wing 
above fuse. C.L. in ft. 
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3. dderiv.m 


% Eric J. Watkiss 

% AA0810 Thesis 

% File to calculate dimensional derivatives 
% dderiv.m 

% Last Update: 12 FEB 94 


% Run nondimensional derivative program 
ndderiv 


% Calculate dynamic pressure 
qbar — .5*rho*Ufps^2; % ft lbs 


Malpha = CMalpha*qgbar*S*cbar/lyy; % per second”2 


Ma = CMg*(cbar/(2*Ufps)*qbar*S*cbar/lyy; 
% per second 


Malphadot = CMalphadot*(cbar/(2*Ufps))*qbar*S*cbar/lyy,; 
% per second 


Xu = -2*CD*gbar*S/(m*Ufps), % per second 
Zu = -2*CL*gbar*S/(m*Ufps); % per second 


Zalphadot = CLalphadot*(cbar/(2* Ufps))*(qbar* S/m); 
% ft per second 


Zq = CLq*(cbar/(2*Ufps))*(qbar* S/m); % ft per second 


Mu = 0; % per ft second 
Xde = -1*CDde*qbar*S/m; % ft per second”2 
Zde = -1*CLdelta*qbar*S/m; % ft per second^2 
Mde = CMde* qbar*S*cbar/Iyy; % per second”2 
Xalpha = (CL - CDalpha)*qbar*S/m, % ft per second^2 


Zalpha = -1*(CLalphaw+CD)*qbar*S/m;  % ft per second”2 
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YB = CyB*qbar*S/m; 

LB = CIB*qbar*S *b/Ixx; 

NB = CnB*qbar*S*b/Izz; 

Yp = Cyp*b*qbar*S/(2*Ufps*m); 
Yr = Cyr*b*qbar*S/(2*Ufps*m); 


Lp 7 Clp*(b/(2*Ufps))*qbar*S*b/Ixx; 


Np 7 Cnp*(b/(2*Ufps))*qbar*S*b/Izz, 


Lr = Clr*(b/(2*Ufps))*qbar*S*b/Ixx; 
Nr = Cnr*(b/(2* Ufps))*qbar* S *b/Izz; 
Y dr = -1*Cydr*qbar*S/m; 

Y da = 0; 

Ldr - Cldr*qbar*S*b/Ixx, 

Lda = Clda*qbar*S*b/Ixx; 

Ndr = Cndr* qbar* S*b/Izz; 


Nda — Cnda*qbar*S *b/Izz; 


% ftísecond”2 


Vo l/second^2 


% l/second”2 


% ft/sec 


% ft/sec 


% 1/sec 


% 1/sec 


% 1/sec 


% 1/sec 


% ft/sec*2 


% ft/sec*2 


9/9 l/sec^2 


% l/sec^2 


9o l/sec^2 


9/9 l/sec^2 


Malphaprime 7 Malpha 4 Malphadot' (Zalpha/Ufps); 


Mqprime — Mq 4 Malphadot; 


Mdeprime = Mde + Malphadot*(Zde/Ufps); 


48 


4. godallas.m 


% Eric J. Watkiss 

9 AA0810 Thesis 

% File to get values of dimensional and nondimensional derivatives 
% godallas.m 

% Last Update: 04 FEB 94 


!del bluebird.dia 
diary bluebird.dia 
diary on 


% Run programs to calculate derivatives 
dderiv 


% Nondimensional Derivatives 
EL 

CD 

CDO 
CLalphaw 
CMalpha 
CDalpha 
CLalphadot 
CMalphadot 
CLq 

CMq 
CLdelta 
CDde 
CMde 

CyB 

CnB 

CIB 


Cyp 
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Cldr 
CDq 
Cyq 


% Longitudinal Dimensional Derivatives 
Xu 

Xalpha 
Xde 

Zu 

Zalpha 
Zalphadot 
Zq 

Zde 

Mu 
Malpha 
Malphadot 
Mq 

Mde 


% Lateral/Directional Dimensional Derivatives 
YB 
Yp 
Yr 
Y da 
Y dr 
LB 
Lp 
Lr 
Lda 
Ldr 
NB 
Np 
Nr 
Nda 
Ndr 


diary off 
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5. ndderiv.m 


% Eric J. Watkiss 

99 AA0810 Thesis 

% File to calculate nondimensional derivatives 
9/5 ndderiv.m 

% Last Update: 12 FEB 94 


% Load Bluebird data with flight condition 
physcalc 


% Calculate coefficients of lift and drag 
CL = W/(.5*rho* Ufps’2*S); 
CD = CL/LD; 


% Calculate lift curve slope of wing in per radian 
CLalphaw = CLalphaafw/(1+CLalphaafw/(pi*ee* AR)); 


% Calculate lift curve slope of horizontal tail in per radian 
CLalphat = CLalphaaft/(1+CLalphaaft/(pi*ee* ARt)); 


% Calculate lift curve slope of vertical tail in per radian 
CLalphav = CLalphaafv/(1+CLalphaafv/(pi*ee* ARv)); 


% Calculate change in hor. tail lift with change in elevator 
dcLtdde = da0dde * CLalphat; % per radian 


% Calculate change in vert. tail lift with change in rudder 
dcLvddr = daOddr * CLalphav; % per radian 


% Calculate zero lift pitching moment 
CMO = CMac + VH * CLalphat * (it + e0); 


% Calculate CMalpha in per radian 
CMalpha = CLalphaw*((h-hac)-VH*(CLalphat/CLalphaw)*(1-deda)); 


% Calculate change in aircraft lift with change in elevator 
CLdelta = dcLtdde*(St/S); % per radian 


% Calculate chng in aircraft pitching moment w. chng in elevator 
CMde --1*VH*dcLtdde; % per radian 


% Calculate angle of attack and elevator angle for trimmed flight 
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% 
% CM = CMO + CMalpha* alpha + CMde*de 
% Cl = CLalphaw*alpha + CLdelta*de 
% 
% EL M ww 
% | CLalphaw CLdelta | | alpha | |CL | 
% | | | a 
% | CMalpha CMde_ || de | | -CMO | 
% 
% A COM NEN 
% 
A = [ CLalphaw CLdelta 
CMalpha CMde ]; 


C=[ CL 
-1*CMO ]; 
X = inv(A)*C; 
atrim = X(1,1); % trim a.0.a. in radians 
etrim = X(2,1); % trim elevator in radians 


% Calculate change in yawing moment with change in rudder 
% "rudder power" 

% assumes VF/Vinfinity = 1 

Cndr = -1* VV *dcLvddr; % in per radian 


% Calculate CnB contribution from vert. tail 

9o CnB 7 CLalphav* VV*(VF/Vinfinity)^2*(1-dsigma/dbeta) 
% assumes VF/Vinfinity = 1 and dsigma/dbeta = 0 

CnB = CLalphav* VV;% in per radian 


% Calculate change in rolling moment with change in sideslip 


% First calculate dihedral contribution from wing 
% Raymer p. 439 

CiBwf = -1.2*sqrt(AR)*Zwf*(Df+W0D/b72; 
CIBw 7 CIBwCL*CL4CIBwf, 


% Next calculate contribution from fin 

yo CIBv 7 -1*Clalphav*Vprime*(VF/Vinfinity)^2*(1-dsigma/dbeta) 
% Assume VF/Vinfinity = 1 and dsigma/dbeta = 0 

CIBv = -1*CLalphav* Vprime; % in per rad 


% Combine ClBg and CIBv into CIB 
CIB = CiIBw + CIBv; % in per rad 
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% Calculate "aileron power", Clda 

9/ See Smetana pp. 139-141 

Cldatau = Cldatauo - Cldataui; 

Clda — Cldatau'tau; % in per radian 


% Calculate change in yawing moment w. aileron deflection 
Cnda - 2*K*CL*Clda; % in per radian 


% Calculate side force due to yaw 
% By Smetana p. 107 
CyB = -.31, % in per radian 


% Calculate side force due to rudder 
Cydr = CLalphav*taur*Sv/S; % in per radian 


% Calculate side force due to aileron 
% By Smetana, p. 138 
Cyda = 0; 


% Calculate rolling moment due to rudder 
Cldr = Cydr*zv/b;  % in per radian 


% Calculate change in drag due to change in elevator 

% Smetana pp. 95-100 

% Using Figure 26 at 8 degrees aoa 

CDde = ((.155-.047)/(20* pi/180))*St/S; % in per radian 


% Calculate change in drag with change in aoa 

“ Smetana pp. 64-65 

9o Assuming dCDO/dalpha is negligible 

CDalpha - 2*CL*CLalphaw/(pi*ee* AR); — Vo in per radian 


% Calculate change in pitching moment w.r.t. alphadot 
% Smetana pp. 78-81, etat assumed = 1 

CMalphadot = -2*CLalphat*deda*(Itprime/cbar)*... 
(It/cbar)*(St/S); % in per radian 


% Calculate change in lift with pitch rate 

% Smetana pp. 82-85 

% Neglecting wing contribution, assuming etat = 1 
CLq = 2*(It/cbar)*CLalphat*(St/S); % in per radian 
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% Calculate change in lift with alphadot 
% Smetana pp. 75-76 
CLalphadot = -1*CMalphadot/(It/cbar); % in per radian 


% Calculate change in pitching moment w. pitch rate 

% Smetana pp. 87-88 

% Assuming etat = | 

CMq = -2*(cbar/4-h)*abs(cbar/4-h)*CLalphaw/(cbar’2) - ... 
2*(It/cbar)^2* CLalphat*(St/S); % in per radian 


% Calculate roll damping 

% Smetana pp. 122-125 

% Neglecting contribution from vertical tail 

Clp = -.475*(AR+4)/(2*pi* AR/CLalphaw+4); % in per radian 


% Calculate change in yawing moment due to rolling 
% Smetana pp. 126-129 

% Neglecting contribution from vertical tail 

Cnp = -1*CL/8; % in per radian 


% Calculate change in side force with yaw rate 
% From Schmidt p. 3-23 

% Assume etat = 1 

Cyr = 2*VV*CLalphav; % in per radian 


% Calculate change in rolling moment w. yaw rate 
% Schmidt p. 3-24 

% Tail contribution 

Clrv = (zv/b)*Cyr; % in per radian 

% Wing contribution 

Clrw = CL/4; % in per radian 

% Total 

Clr » Clrv * Clrw;  % in per radian 


% Calculate yaw damping 

% Schmidt p. 3-25 

% Tail contribution 

Cnrv = -1*(lv/b)*Cyr; % in per radian 

% Wing contribution from Smetana p. 136 
CDO = CD-CL“2/(pi*ee* AR); 

Cnrw = -.02*CL‘%2-.3*CD0; % in per radian 
% Total 

Cnr 7 Cnrv * Cnrw; % in per radian 
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% The following 3 derivatives are negligible and taken to be O 
CDq = 0; % in per radian 
Cyq=0; % in per radian 
eyp = 0; % in per radian 


% A few misc. calculations 
% Static Margin/Neutral Point 


statmar = CMalpha/(-1*CLalphaw) 
hn = statmar + h 
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6. physcalc.m 


% Eric Watkiss 

Vo AA0810 Thesis 

% File to calculate physical considerations 
% physcalc.m 

% Last Update: 04 FEB 94 


% Load fixed Bluebird data 
bluebird 


% Load flight condition 
bibrdfc l 


% Calculate aircraft weight 
W = WImg + Wrmg + Wng; 


% Calculate aircraft mass 
m = W/g; 


% Calculate aspect ratio of wing 
AR 7 b^2/S; 


% Calculate aspect ratio of hor. tail 
ARt — bt^2/St, 


% Calculate aspect ratio of vert. tail 
ARv = bv“2/Sv; 


% Calculate longitudinal center of gravity 
h = ((ng*Wng + mg*(Wlmg+Wrmg)/W-lewing)/cbar, 


% Calculate "tail length" from c.g. to horizontal tail a.c. 
% same for horizontal and vertical 

It = c4tail - (lewing + h*cbar); 

lv = It; 


% Calculate "tail length" from c/4 wing to c/4 tail 
Itprime = c4tail - c4wing; 


% Calculate hor. tail volume coefficient 
VH = It*St/(S*cbar); 
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% Calculate vert. tail volume coefficient (yaw) 
VV = lv*Sv/(b*S); 


% Calculate vert. tail volume coefficient (roll) 
Vprime = zv*Sv/(b*S); 


% Unit antisymmetrical angle of attack for outer and inner 
% edge of aileron (See Smetana p. 141) 

antisymo = ao/(b/2); 

Cldatauo = .74; 

antisvmi — al/(b/2); 

Cldataui = .23; 

cacw = ac/cbar; 

tau = .52; 


% for yawing moment due to aileron, see p. 142, Smetana 
eta — ai/(b/2); 
K = -.17; 


% for side force due to rudder deflection, see Smetana p. 145 
vratio = Sr/Sv; 


taur = .62; 


Vo for rolling moment due to sideslip, See Raymer, Fig. 16.21, p. 439 
CIBWCL = -.04; 
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APPENDIX B: MOMENT OF INERTIA CALCULATIONS 


For the formula 


5 52 
_ [Wa + 20(L short + LLoNG)Ilm:5 Fu Z M 
AA y 


4n? $ 


20| r2 Z 
SS (L SHORT + 2 


the following data were used for extraction of the moments of inertia. 


Fixed values: 


Weight of the model, 

Weight per unit length of chain, 
Length of short chain, 

Length of long chain, 
Gravitational constant, 


(adjusted for latitude and elevation) 


Variable values: 


ix 


Distance from pivot to center of gravitv 
of model and support, 

Distance from pivot to center of gravity 
of model, 

Average period of model and support, 


Distance from pivot to center of gravity 
of model and support, 

Distance from pivot to center of gravity 
of model, 

Average period of model and support, 


Distance from pivot to center of gravity 
of model and support, 

Distance from pivot to center of gravity 
of model, 

Average period of model and support, 
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Wu = 58.45 lbs 

G — 0.06133 Ibs/ft 
Lshorr = 13.00 ft 
Lione 2 15.00% 
g = 32.1472 ft/sec’ 


Zuss = 12.95 ft 


Lu USA 
Puss = 4.109 sec 


12.01 ft 


Zug 
Zu = 12358 
Pu+s = 3.976 sec 


Zas = 12.49 ft 


Zu = 12.81 ft 
Puss = 4.077 sec 


APPENDIX C: MATLAB PROGRAMS USED TO ESTIMATE 
AIRCRAFT DYNAMICS 


A. LONGITUDINAL DYNAMICS 
l.  longnat.m 


% Determines the short and long period natural response 
% using the full longitudinal plant 


clear 

dderiv 

Idel longnat.dia 
diary longnat.dia 


% Inertial Matrix 

in = [Ufps 0 0 0; 
O  Ufps-Zalphadot 0 0; 
0 -1*Malphadot 1 0; 
0 0 O 1); 


90 Aircraft Matrix 
an —[Ufps*Xu Xalpha 0 -1*g*cos(thetanaut); 
Ufps*Zu Zalpha Ufps+Zq  -1*g*sin(thetanaut); 
Ufps*Mu Malpha Mq 0; 
0 0 l 0]; 


a — inv(in)'an; 

b =(1,1,1,1) 

c =eye(size(a)); 

d = [0; 0; 0; 0], 

t = 0:.01:100; 

[y,x,t] = impulse(a,b,c,d, 1,t); 
clg 


figure(1) 
Yoplotting alpha and q, short period 
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plot(t,y(:,2));axis([0 15 -1.5 1.5));grid;gtext('alpha') 
hold on 

plot(t,v(:,3));gtext('q') 

xlabel('Time, sec');vlabel(' Response Gain') 
%title(‘Natural Response’) 

hold off 

pause 


figure(2) 

% plotting u/U and theta, phugoid 

plot(t,y(:, 1));axis({O 100 -2 2]);grid;gtext(‘u/U’) 
hold on 

plot(t,y(.,4));gtext(‘theta’) 

xlabel('Time, sec');ylabel('Response Gain”) 
%title( "Natural Response”); 

hold off 


damp(eig(a)) 


diary off 
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Z. stepper.m 


% Determines the short and long period step response 
% using the full longitudinal plant 


clear 
dderiv 


% Inertial Matrix 

in — [Ufps 0 0 0. 
O Ufps-Zalphadot O 0; 
O -1*Malphadot 1 0; 
0 0 0 1; 


% Aircraft Matrix 

an —[Ufps*Xu Xalpha O -l*g*cos(thetanaut); 
Ufps*Zu Zalpha Ufps-Zq -1*g*sin(thetanaut); 
Ufps*Mu Malpha Maq O; 
0 0 l 0); 


% Control Matrix 
bn-[Xde Zde Mde 0]' 


a — inv(in)'an; 
b — inv(in)'bn; 
c — eve(size(a)); 
d -[0000]; 

t = 0:.01:100; 


[y,x,t] = step(a,b,c,d, 1,t); 


clg 

figure(1) 

Yplotting alpha and q, short period 
plot(t,y(:,2)) 

axis({0 15 -5 3]) 

grid;gtext('alpha)) 

hold on 

plot(t,y(:,3));gtext('q') 

xlabel('Time, sec');vlabel(Response Gain") 
%title('Short Period Step Response") 
hold off 
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pause 


figure(2) 

%plotting u/U and theta, phugoid 
plot(t.yC,1)) 

axis 

grid;gtext('u/U') 

hold on 

plot(t, y(:,4));gtext(theta') 

xlabel( Time, sec);ylabel('Response Gain') 
%title('Long Period Step Response"); 
hold off 

pause 


figure(3) 

Veplotting alpha and q, long period 
plot(t,y(-,2)) 

grid;gtext('alpha') 

hold on 

plot(t,v(,3));gtext('q') 

xlabel(' Time, sec');vlabel('Response Gain') 
%title('Long Period Step Response”) 

hold off 
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3. n step.m 


% This program uses STEP function to determine 
% normal acceleration response to a unit step input 


clear 
dderiv 


% Inertial Matrix 
= [Ufps 0 0 0; 
O Ufps-Zalphadot O 0; 
O -1*Malphadot 1 0; 
0 0 0 1]: 


% Aircraft Matrix 
=[Ufps*Xu Xalpha 0 -1*g*cos(thetanaut); 
Ufps*Zu Zalpha Ufps+Zq -1*g*sin(thetanaut), 
Ufps*Mu Malpha Mq 0; 
0 0 1 0]; 


% Control Matrix 
bn = [Xde Zde Mde 0]'; 


a =inv(in)*an; 

b = inv(in)*bn; 

c =[ Zu Zalpha Zq 0}; 
d =Zde; 

t =0:0.05:15; 


[y.x,t] = step(a,b,c,d, 1,t); 


for 1 = 1:length(t) 
y2 = y/32.174* pi/180; 
end 


% plotting g's/deg 

plot(t,y2);grid 

xlabel('Time, sec');ylabel('Normal Acceleration: g s/deg') 
%title('Acceleration Response to Unit Step') 
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4.  homogen.m 


% Determines the short and long period homogeneous response 
% using the full longitudinal plant 


clear 
dderiv 


% Inertial Matrix 

in = [Ufps 0 0 0 
O Ufps-Zalphadot 0 
O -1*Malphadot 1! 
0 0 o 1); 


> 


0; 
Q: 


> 


% Aircraft Matrix 

an =[Ufps*Xu Xalpha O -1*g*cos(thetanaut); 
Ufps*Zu Zalpha Ufps+Zq -1*g*sin(thetanaut); 
Ufps*Mu Malpha Mq 0; 
0 0 l 0]; 


% Control Matnx 
bn =[Xde Zde Mde O]’; 


a =inv(in)*an; 


b =inv(in)*bn; 
c =eye(size(a)); 
d -[0000] 

t = 0:.01:100; 


[y,x,t] = step(a,b,c,d,1,t); 


x0 = y(find(t==15),:); 
XU 7 x0/x0(2)* 5/57.3 
[y,x,t] = initial(a,b,c,d,x0,t); 


clg 

figure(1) 

%plotting alpha and theta, short period 
plot(t,y(:,2)) 

%title('Short Pernod Homogeneous Response’); 
grid 

axis([0 15 -.2 .2)); 
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gtext('alpha - rad'); 

hold on 

plot(t,y(:,3)),gtext(q - rad/sec'); 
xlabel('Time, sec”); ylabel('Response Gain”), 
hold off 

pause 


figure(2) 

Y%plotting u/U and theta, long period 
plot(t,y(:,1)) 

Votitle(Long Period Homogeneous Response’); 
grid 

axis; 

gtext('u/U'); 

hold on 

plot(t,y(:,4));gtext(theta - rad"); 
xlabel('Time, sec'), vlabel( Response Gain'); 
hold off 
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5. doublet. m 
% This program uses EXPM function to find response to pitch doublet 


clear 
dderiv 


% Inertial Matrix 
= [Ufps 0 O O; 
O Ufps-Zalphadot 0 0; 
O -1*Malphadot 1 0; 
0 0 o 1); 


% Aircraft Matrix 
=[Ufps*Xu Xalpha O -1*g*cos(thetanaut); 
Ufps*Zu Zalpha Ufps+Zq -1*g*sin(thetanaut); 
Ufps*Mu Malpha Mq 0; 
0 0 l 0]; 


bn = [ Xde Zde Mde 0 ]'; 


a =inv(in)*an; 
b -inv(in)*bn; 


tl =1.0; % start of doublet, sec 
t2 =2.0; % midpoint of doublet, sec 
t3 23.0; %end of doublet, sec 


dO —-1; Vo unit elevat. nput (1 rad) [t.e.u] 
tim 7 15; % set end time 
t =0:0.05:tim; 


x = zeros(4,length(t));  % initialize solution matrices 
x1 = zeros(4,length(t)); 
x2 = zeros(4,length(t)); 
x3 = zeros(4,length(t)); 


for 1 = 1:length(t) 


if t(1) >= tl 
de(1) - dO; 
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x1(,1) = d0*(expm(a*(t(1)-t1)) - eye(size(a))*inv(a)*b,; 
end 


if t(i) >= t2 

de(i) = de(1) - 2*d0; 

x2(:,i) = -2*d0*(expm(a*(t(i)-t2)) - eye(size(a)))*inv(a)*b; 
end 


if t(1) >= t3 

de(i) = de(1) + dO; 

x3(:,1) = dO*(expm(a*(t(1)-t3)) - eye(size(a)))*inv(a)*b; 
end 


x(,1) = x1(:,1) * x2(,1) * x3(:,1); 
end 


ymin — -2*abs(dO); 
ymax = 2*abs(d0); 
V=[0 tim ymin ymax]’; 


figure(1) 

plot(t,de);grid;axis(V) 

xlabel('Time, sec'); ylabel('Elevator, Rad’); 
Vetitle('Elevator Input") 

pause; axis 


figure(2) 

Voplotting alpha and pitch rate 
plot(t,x(2,:)) 

“%title(‘Doublet Response’) 
grid;gtext('alpha') 

hold on 

plot(t,x(3,:));gtext('q') 

xlabel('Time, sec');vlabel('Response Gain') 
hold off 
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6. sp bode.m 


% This program uses the BODE function to find the 
%  short-period response (transfer-function gain and phase) 
% to a harmonic input 


clear 
dderiv 


% Inertial Matrix 
= [Ufps 0 0 0; 
O Ufps-Zalphadot O 0; 
O -1*Malphadot 1 0; 
0 0 OEM: 


% Aircraft Matrix 
= [Ufps*Xu  Xalpha 0 -1*g*cos(thetanaut); 
Ufps*Zu Zalpha Ufps+Zq -1*g*sin(thetanaut), 
Ufps*Mu Malpha Mq 0; 
0 (e 0): 
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bn-[Xde Zde Mde 0]; 


a =inv(in)*an; 

b = inv(in)*bn; 

c =eye(size(a)), 
d — zeros(size(b)); 


= logspace(-1,2,150); 
[mag, phase, w] = bode(a,b,c,d, I, w); 


for 1 = 1:length(w) 


for j= 1:4 
if phase(i,j) > 0 
phase(i,j) = phase(1,j) - 360; 
end 
end 
end 
clg 
figure(1) 


% Plotting alpha and pitch-rate gains 
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semilogx(w,mag(:,1));grid;gtext('alpha') 
%title('Short-Period Frequency Response’) 
hold on 

semilogx(w,mag(:,2)),gtext('q') 
xlabel('Frequency, rad/sec');ylabel('Gain') 
hold off 

pause 


figure(2) 

% Plotting alpha and pitch-rate phase angles 
semilogx(w,phase(:,1));grid;gtext('alpha") 
%title('Short-Period Frequency Response”) 
hold on 

semilogx(w,phase(:,2));gtext('q') 
xlabel('Frequency, rad/sec');vlabel('Phase, deg 
hold off 
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7. n bode.m 


% This program uses the BODE function to find 
% the normal-acceleration response 
%  (transfer-function gain and phase) to a harmonic input 


clear 
dderiv 


% Inertial Matrix 

in = [Ufps 0 OO: 
O Ufps-Zalphadot O 0; 
O -1*Malphadot 1 0; 
0 0 0,3 


90 Aircraft Matrix 
— (UfpstXu  Xalpha 0 -1*g*cos(thetanaut); 
Ufps*Zu Zalpha Ufps+Zq -1*g*sin(thetanaut); 
Ufps*Mu Malpha Mq O, 
0 0 l 0]; 


% Control Matrix 
bn = [Xde Zde Mde 0)'; 


a =inv(in)*an; 

b -inv(in)*bn; 

c =[ Zu Zalpha Zq 0 ]; 
d =Zde; 


= logspace(-1,2,150); 
[mag,phase,w] = bode(a,b,c,d, 1,w); 


for i = 1:length(w) 
mag(i) = mag(i)/g*pi/180; % converting to g's/deg 
if phase(i) > 0 
phase(i) = phase(1) - 360; 
end 
end 


% Plotting load factor/deg elevator input 

figure(1) 

semilogx(w,mag);grid 

xlabel('Frequency, rad/sec'); ylabel("Normal Acceleration Gain’) 
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%title(‘Normal Acceleration Response to Harmonic Input’) 
pause 


% Plotting normal acceleration phase angle 

figure(2) 

semilogx(w, phase); grid 

xlabel('Frequency, rad/sec'); ylabel('Normal Acceleration Phase, deg’) 
Voetitle "Normal Acceleration Response to Harmonic Input) 
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B. LATERAL-DIRECTIONAL DYNAMICS 
l. lat dir.m 


% Solves the full 4x4 lateral-directional response 
clear 
dderiv 
!del latdir.dia 
diary latdir.dia 
in -[Ufps 0 0 0 
0 l O -1*Ixz/lxx 
UO l 0 
0 -I*Ixz/Izz O 1) 
an-[YB Yp g*cos(thetanaut) Yr-Ufps 


LEER 0 Lr 
0O 1 0 0 
NB Np 0 Nr]; 


a = inv(in)*an; 
[wn,zeta] = damp(a) 
[x,d] = eig(a) 

r=(0 0 0 OJ; 


fori= 1:4 
U) = d(1,1); 


end 
x2 = zeros(4,2); 


forj=1:4 

x2(], 1) 7 abs(x(j,1)); 

x2(j,2) = 180/pi*angle(xG, 1)); 
end 


xx = x2(1,1); 

xy = x(3,2); 

xz = x(4,3); 

x3(:,1) = x2(:,1)/xx; 
x3(:,2) = x2(:,2) - x2(1,2); 
x3(:,3) = x(,3)/xy; 
X3(:,4) — x(:,4)/xz; 

x3 


diary off 
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2.  rudkick.m 
% Response to unit rudder impulse 


clear 
dderiv 


in-[Ufps 0 0 0 
B 1 O -1*Ixz/Ixx 
0 0 l 0 
O -1*lxz/Izz 0 1]; 


y/o Plant Matrix 

an=[YB Yp g*cos(thetanaut) Yr-Ufps 
LB Lp 0 Er 
O 1l 0 0 
NB Np 0 Nr]; 


% Control Matrix 
bn = [Ydr Yda 
Ldr Lda 
0 0 
Ndr Nda], 


a — inv(in)*an; 

b = inv(in)*an; 

c = eye(size(a)); 

d = zeros(size(b)); 


t 70:05:10; 
[x,y,t] = impulse(a,b,c,d, 1,t); 


% plotting sideslip (beta) and bank angle (phi) 
plot(t,v(:,1));grid;gtext('beta') 

hold on 

plot(t,v(,3));gtext('phi') 

xlabel('Time, sec');vlabel(' Response Gain!) 
Votitle('Rudder Kick Response') 

hold off 

pause 


% integrate p for final bank angle 
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phi = 0; 
t =O 05100; 


[x,y,t] = impulse(a,b,c,d, 1,t); 
for 1 = l-length(t) 
phi = phi + y(i,2)*.05; 


end 


phi 
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3.  sideslip.m 


% Solves the full 4x4 lateral-directional response for 
% initial sideslip condition 


clear 
dderiv 


in=[Ufps 0 0 0 
| O -I*Ixz/Ixx 
"m 0 l 0 
O -1*Ixz/Izz 0 1]; 


an=[YB Yp g*cos(thetanaut) Yr-Ufps 
LB Lp 0 Lr 
n ] 0 0 
NB Np 0 Nr]; 


bn-[Ydr Yda 
Ldr Lda 
DE 0 
Ndr Nda]; 


a — inv(in)*an; 

b = inv(in)*bn; 

c = eye(size(a)); 

d — zeros(size(b)); 


t = 0:.05:15; 
% calculate sideslip per bank angle 


f-(00-IĊCL f' 
k=[ CnB Cndr Cnda 
CIB Cldr Clda 

CyB Cydr Cyda J; 


perphi = inv(k)*f; 
betaperphi = perphi(1); 


PHI - 10*pi/180; 


BETA - betaperphi* PHI, 
xin - [ BETA 0 PHI 0J' 
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[x,y,t] = initial(a,b,c,d,xin,t); 


plot(t,y(:, 1), -);grid;gtext('beta") 

%title(Dutch Roll Response for Sideslip I. C.s') 
hold on 

plot(t,y(:,3),'~');gtext(‘phi’) 

hold off 

xlabel('Time in Seconds”) 

ylabel(' Amplitude”) 
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4. roll. m 


% Solves the full 4x4 lateral-directional for roll response 
% (transfer function gain and phase) due to 
% a harmonic aileron input 


clear 
dderiv 


in= [Ufps 0 0 0 
no] O0 -1*Ixz/Ixx 
nm 0 l 0 
0 -I*Ixz/Izz O 1] 


an-[YB Yp g*cos(thetanaut) Yr-Ufps 
LB Lp 0 Lr 
0 1 0 0 
NB Np 0 Nr] 


bn-[Ydr Yda 
Ldr Lda 
0O 0 
Ndr Nda] 


a = inv(in)*an 

b = inv(in)*bn 

c — eye(size(a)) 

d = zeros(size(b)) 


w = logspace(-1,2,150); 
[mag,phase,w] = bode(a,b,c,d,2,w); 


for 1 = 1:length(w) 
for j= 1:4 
if phase(i,j) 9 0 
phase(i,j) — phase(1,]) - 360; 
end 
end 
end 


clg 


figure(1) 
% Plotting roll angle gain 
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semilogx(w,mag(:,3));grid 

Votitle'Harmonic Aileron Input") 

xlabel( Frequency, rad/sec^);ylabel('Response Gain') 
pause 


figure(2) 

% Plotting roll angle phase 
semilogx(w,phase(:,3));grid 

Votitle'Harmonic Aileron Input') 

xlabel( Frequency, rad/sec');vlabel('Phase in Degrees') 
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